function J_i = HJB_t_OJS( J_old, param, aux)

    ind_stay = 1   + (aux.n_m+1).*(0:aux.n_m-1);

    T_tr=param.T+param.T_lam.*(param.s.*param.lambda+param.zeta);
    A=1./aux.dt*speye(aux.n_m,aux.n_m)+aux.CN.*T_tr;
    B=1/aux.dt.*speye(aux.n_m,aux.n_m)-(1-aux.CN).*T_tr;

    J_i=J_old;
    J_i_old=J_i;

    diff_J=1;
    u=1;
    while diff_J>=10^(-16)

        income   = (1-param.omega_1).*exp(param.ln_m_HJB) -param.omega_0 + B*J_old;
        P_m=aux.pen.*(J_i<0);
        P_p=aux.pen.*(J_i>param.c);

        AA=A;
        AA(ind_stay)=AA(ind_stay)'+P_m+P_p;
        J_i=AA\(income+P_p.*(param.c)+P_m.*(0));


        diff_J=max(abs(J_i_old-J_i));
        J_i_old=J_i;
        u=1+u;
        if u>10^4
           error('HJB failed to coverge')
        end
    end
end

